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ABSTRACT 

The influence of rotation on the dynamical evolution of collisional stellar systems is 
investigated by solving the orbit-averaged Fokker-Planck equation in {E, Jz)-s\ia.ce. 
We find that large amounts of initial rotation drive the system into a phase of strong 
mass loss while it is moderately contracting. The core is rotating even faster than 
before although angular momentum is transported outwards. At the same time the 
core is heating. Given these features this phase can be associated with the gravo-gyro 
'catastrophe' found by Hachisu (1979). The increase in central angular momentum 
levels off after about 2-3 initial half-mass relaxation times indicating that the source 
of this 'catastrophe' is depleted. Finally, the central angular velocity increases again, 
but with a rather small power of the central density - the same power as for the 
central velocity dispersion during self-similar contraction towards the collapse. The 
rotation curves flatten and the ellipticity variations decrease with time, but their 
shapes are very similar. These results suggest the existence of a self-similar solution 
for a rotating cluster as well. The maximum values of rotational velocity and ellipticity 
occur at about the half mass radius. 

Key words: methods: numerical - celestial mechanics, stellar dynamics - globular 
clusters: general. 



1 INTRODUCTION 

Tlie study of tlie dynamical behaviour of spherically sym- 
metric, collisional stellar systems has made considerable 
progress during the last decade. Numerical improvements 
and refinements are applied to current techniques (Fokker- 
Planck, Monte-Carlo, fluid-dynamical or N-Body) in order 
to investigate anisotropy and the post-collapse phase (Bett- 
wieser & Sugimoto 1984, Murphy et al. 1990, Heggie & 
Ramamani 1989, Takahashi 1996a,b, Spurzem 1994, 1996), 
stochastic energy input due to binaries (Giersz 1996) , small- 
N statistics (Giersz & Heggie 1994ab, Heggie 1996), primor- 
dial binaries (Hut 1996), tidal effects, stellar evolution or 
mass spectra (Chernoff & Weinberg 1990, Lee et al. 1991, 
Fukushige & Heggie 1995, Drukier 1995, Giersz & Heggie 
1996, 1997). Considering anisotropy, there are still other fea- 
tures in addition to its hitherto studied effects on the radial 
variations of density and velocity dispersion, one important 
effect being the feasibility to flatten those systems. Binney 
(1978) applied the tensor virial theorem to elliptical galax- 
ies constructed from similar ellipsoids and comprising flat 
rotation curves, and he derived formulae relating ellipticity 
to the ratio of rotation velocity and velocity dispersion. The 
observed spread in Vm/cr for galaxies of a given ellipticity 



is explained as an effect of anisotropy, which changes the 
hydrostatic equilibrium. 

Observations show that flattening is a common feature 
of globular clusters, which has been known since the early 
work done by Pease & Shapley (1917). Measuring projected 
ellipticities e — 1 — b/a of large globular cluster samples 
White & Shawl (1987) derive a mean e = 0.07 ± 0.01 for 
99 clusters in the Milky Way, and Staneva et al. (1996) flnd 
e — 0.086 ± 0.038 for 173 clusters in M31, with maximum 
values 0.27 and 0.24 of individual globulars, respectively. 
Kinematical data, i.e. radial velocities of large numbers of 
cluster members, reveal that this flattening may indeed be 
explained in terms of rotation, and that the minor axes are 
nearly coincident with the determined rotation axes (Mey- 
lan & Mayor 1986). Dust obscuration, anisotropy or tidal 
distortion are able to explain individual cases of flattening, 
but can statistically be ruled out as the main mechanism 
(White & Shawl 1987). 

Significant ellipticity variations are found within glob- 
ular clusters (e.g., Geyer & Richtler 1981, Geyer et al. 
1983), and these partly also coincide with the rotation 
curves obtained by fits to the radial velocity data with some 
parametrization specifled for the velocity field (Meylan & 
Mayor 1986). Moreover, Kontizas et al. (1990) show that 
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the outer parts of globulars in the Small Magellanic Cloud 
are obviously rounder than the parts inside the half mass 
radius and it is likely that their structure differs from that 
of the galactic globular clusters because they are younger 
(in general) and subject to different tidal forces. The impor- 
tance of age for the interpretation of observed ellipticities 
has already been emphasized by Frenk & Fall (1982), who 
undertook eye-estimates of cluster ellipticities in the Milky 
Way and the Magellanic Clouds, the latter being slightly 
larger than the former, which is explained again in terms of 
internal globular cluster evolution. This view is supported 
by studies relating Milky Way globular cluster ellipticities to 
the cluster concentration parameter c — log{rt/rc) (White 
& Shawl 1987), where rt is the tidal radius and Vc is the 
core radius, or to the half mass relaxation time (Davoust 
& Prugniel 1990), both representing the dynamical age of a 
globular cluster. In these two investigations, the average flat- 
tening of the dynamically younger systems is significantly 
larger as well, indicating that loss of angular momentum, 
presumably originating from diffusion past the escape ve- 
locity on relaxation time scales, decreases the ellipticity of 
a cluster. 

Indeed, Agekian (1958) suggested a model in which spe- 
cific angular momentum is lost due to a relatively large 
fraction of escaping stars residing in the tail of a rotat- 
ing Maxwellian velocity distribution shifted towards the di- 
rection of rotation as compared to the fraction which is 
counter-rotating. Considering this effect for every volume 
element of rotating ellipsoids he obtained a critical elliptic- 
ity of e = 0.74 below which the systems become rounder 
with time. 

More recently, even larger compilations of radial veloc- 
ity data have been obtained by different working groups 
for some globulars. With these data sets they are able to 
model the velocity field around the core with nonparamet- 
ric methods, and find in the core of 47 Tuc increasing an- 
gular velocities towards the centre (Gebhardt et al. 1995), 
which is interpreted as a gradually faster rotating core due 
to (gravothermal) contraction of an initially slowly rotat- 
ing solid body. This would imply that trapping of angular 
momentum in the contracting core wins over angular mo- 
mentum transport due to viscosity effects. 

Analogously to the discussion whether evaporation or 
gravothermal contraction drives cluster evolution predomi- 
nantly, a new question arose due to the stability analysis ap- 
plied to adiabatically confined, rotating cylinders by Hachisu 
(1979). He found that, depending on central temperature 
and angular velocity, there exist cases where, if angular mo- 
mentum is removed from a shell, gravitational contraction 
results in an increase in angular velocity, which leads to a 
runaway in angular momentum transport (by viscosity) and 
central contraction. In analogy to 'gravothermal catastro- 
phe' this effect is called 'gravo-gyro catastrophe', and the 
physical origin in the latter case is due to a negative specific 
moment of inertia, analogously to negative specific heat in 
the former case. In order to study these effects accurately 
in a cluster- like structure Akiyama & Sugimoto (1989) set 
up N-Body simulations for rotating clusters, however, they 
simulated 1000 particles only. Therefore, the statistical qual- 
ity of their data is rather poor and consequently, the four- 
phase evolution proposed by them consisting oi i) violent re- 
laxation, a) gravo-gyro instability. Hi) static evolution, iv) 



gravothermal collapse, seems to be rather suggestive than 
indicative. 

Goodman (1983) performed the first numerical simula- 
tions solving the orbit-averaged Fokker-Planck equation in 
(-E, Jz)-space applied to rotating stellar systems. His main 
results are that the difference in ellipticities between old 
galactic globular clusters and the younger Magellanic clus- 
ters can be explained by relaxation processes, though runs 
representing appropriate cluster structure parameters (e.g., 
smaller King-Wo-parameters than 9.8) were not stated in 
his thesis, and that during the self-similar collapse phase the 
central angular velocity increases as a small power of cen- 
tral density, which is explained by linear perturbation anal- 
ysis performed on the self-similar solution found by Lynden- 
Bell and Eggleton (1980). However, this work suffered from 
coarse numerical resolution and from the problem that de- 
terminations of asymptotic values remained open. 

A sequence of rotating King models is established by 
Longaretti & Lagoute (1996), who present an evolutionary 
path represented by the time evolution of structural param- 
eters; their models take into account evaporation and loss of 
angular momentum by using simplified expressions and ap- 
proximations for the derivation of their diffusion coefficients. 
Some of them emerge to be strong constraints to the general 
view of evolutionary characteristics of coUisional sytems. In 
particular, deviations from a rotating King form of the dis- 
tribution function may be significant in advanced stages of 
cluster evolution. 

The present study is aimed to disentangle the effects of 
evaporation and angular momentum transport stated above 
and investigating the detailed evolutionary scenarios, when 
initial models are given, thereby continuing the work begun 
by Goodman (1983). Moreover, the incorporation of the dy- 
namical effects of rotation is seen as one of the stepping 
stones towards realistic models of globular clusters; further 
steps as they have been studied in the non-rotating case 
(post-collapse, mass spectrum, stellar evolution) are the sub- 
ject of future papers. A new self-consistent Fokker-Planck 
solver has been written and is presented together with de- 
tailed numerical results. Section ^ gives a description of how 
the problem is modelled by the Fokker-Planck equation in an 
axisymmetric coordinate system together with the Poisson 
equation, while a detailed derivation of the diffusion coef- 
ficients is given in the Appendix. Also, possible limitations 
of the method due to a neglect of the third integral will be 
discussed. Section ^ provides an overview of the numerical 
method, and in Section ^ numerical results of the present 
single mass models are presented. Section |^ summarizes the 
derived effects of rotation on dynamical evolution, relates it 
to work from other authors and gives prospects for future 
work. 



2 ORBIT AVERAGED FOKKER-PLANCK 
EQUATION 

2.1 Method 

The method described here may be regarded as a continu- 
ation of Goodman's (1983) work and the mathematical for- 
mulation used here is very near to his, but for the reason 
that it belongs to a part ("Paper III") of his Ph.D. -thesis. 
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which was never published elsewhere, a description will be 
repeated here. 

In the present study we follow the evolution of the dis- 
tribution function / as a function of the energy E and the z- 
component of the angular momentum Jz, both representing 
velocity variables, with time t. We consider E and Jz as the 
only isolating integrals and neglect evident non-ergodicity 
on the hypersurface in phase space given by E and Jz due 
to any third integral. Because the derivation of the diffu- 
sion coefficients necessitates an isolating integral expressed 
as a function of coordinate space and velocity variables, an- 
alytical expansions turn out to be unsuitable for the present 
technique and as an approximate third integral in the lim- 
iting case of spherical symmetry would be the only choice 
left to be considered, as it, e.g., has been used in the approx- 
imate self consistent dynamical models prepared by Lupton 
& Gunn (1987). But at the current status of numerical re- 
sources we decided to leave the investigations of 3-integral 
models for the future. 

The Boltzmann equation transformed to the new vari- 
ables then reads as 



21 d^dl 
dt dt dE 



(1) 



with the potential (p advanced according to the Poisson equa- 
tion 

= iTvGn , (2) 

and the coUisional term on the right hand side of Eq. ex- 
pressed under the Fokker-Planck assumption of small angle 
scatterings 



(dl) 



^^«AE>fV) 
-^(<AJ.>m 



(3) 



+ 



-(< ASA Jz > fV) 



dEdJ, 



where V is the volume element given by 2tt/p, with p be- 
ing the radius in cylindrical coordinates. Because the relax- 
ation time tr is much longer than the dynamical time tdyn 
of the system, the Boltzmann equation in our problem re- 
duces to an equation in which the distribution function only 
depends on orbital constants and time. Thus we are allowed 
to take the orbit average over that area in the meridional 
plane which intersects with the hypersurface in phase space 
for which E and Jz are specified. The condition is that ki- 
netic energy in the meridional plane be non-negative: 



The volume of this hypersurface is given by 



p{E, Jz) = 4n^ 



dpdz . 



(4) 



(5) 



A{E,J^) 



where the intersection defined by Eq. 1^ is given by 
A{E,Jz), and the factor in front of the integral is due 
to integration over the third velocity variable, e.g. = 



arctan(up/iiz). Symmetry about the azimuthal direction in 
coordinate space was assumed. 

The orbit averaged Boltzmann equation then takes the 
form 



a/ lag 9/ 

dt pdt dE 



in which we introduce the function q{E, Jz) defined to be 



g(S, Jz) = 2tv' 



(vp + fz) dpdz ; 



(7) 



A(E,J^) 



note that the integrand is equal to 2{E — (/;) — {J^ / p^) and 
evaluates to zero on the boundary of A. Goodman empha- 
sizes that q is the axially symmetric analogue of the radial 
action 



Q{E,J)^2 



Vrdr 



(8) 



= 0, 



(9) 



in a spherically symmetric potential and comprises actually 
the average of the radial action over J with boundaries Jz 
and Jmax, if the potential is spherically symmetric, thus rep- 
resenting an average of an adiabatic invariant with respect 
to the third integral. 

Interpreting q{E,Jz) as an adiabatic invariant itself, 
equation ^ indicates that, neglecting any encounters, 

from which it follows, that there is a redistribution of ener- 
gies in the system, but the 'adiabatic invariant' q and the an- 
gular momentum are conserved. The desired solution of the 
Boltzmann equation therefore may be split into two parts: 
the first is to calculate the evolution due to stellar collisions 
(i.e. small angle scatterings: the Fokker-Planck step) with 
the potential held fixed and the second is to advance the 
distribution function / due to slow adiabatic changes in the 
potential with f{q, Jz) = const. 

The Fokker-Planck equation is transformed into flux 
conservation form in order to improve conservation of sev- 
eral quantities: 

df _1 ( OFe _ dFj^ 
^~dE 



dt 



P 



dJz 



(10) 



with particle flux components 
Fe 

df 



DEEg^-DEJ.gJ- 



DEf 



F. 



(11) 
(12) 



The orbit averaged flux coefficients Da in these equations 
are derived from the local diffusion coefficients, which again 
are found for the axial symmetric geometry via the prescrip- 
tions of Rosenbluth et al. (1957) involving covariant deriva- 
tives instead of the procedure employed by Goodman (1983) 
using a non-covariant form. The derivation of the flux coef- 
flcients is given in the Appendix. Comparison of the expres- 
sions for the diffusion coefficients found by us and Goodman 
reveals complete agreement. 

The derivation of the diffusion coefficients makes it 
necessary to specify a background distribution function by 
which test stars are scattered. Using self-consistently the 
foreground distribution for the background would require 
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computational time proportional to Np x Nz x N% x Nj^ , 
where the Ni are grid sizes in coordinate and velocity space, 
which would reduce the grids employed to inaccurate de- 
scriptions of the current problem; thus, as in all previous 
applications concerning 2D Fokker-Planck methods for stel- 
lar systems, we set up an appropriate form of the background 
distribution. Herein we follow Goodman (1983) giving a ro- 
tating Maxwellian velocity distribution to the background, 



3- exp(- 



2cr2 



) 



(13) 



where p, Q and a correspond to the zeroth, first and second 
order moments of the distribution function, and are the den- 
sity, angular velocity and one-dimensional velocity disper- 
sion of the field star distribution, respectively. Applying this 
form (Eq. (^^) to the background distribution yields ana- 
lytical expressions in the diffusion coefficients parametrized 
locally by the three moments mentioned above, just thereby 
reducing the computational efforts necessary. 



2.2 Initial conditions 

For starting configurations rotating King models are uti- 
lized. These are of the form 

U{E, .h) = const ■ (e-"^ - 1) • e--^""^- , (14) 

where /9 = l/(m(j^) and the dimensionless angular veloc- 
ity LOo = ^J%lAi^Gpc, ■ ilo are parameters to be specified for 
each model. In order to construct a potential-density pair, 
jj is related to the King-parameter (dimensionless potential) 
Wo = — f3m{<l)— <l)t) ■ Throughout this paper our different evo- 
lutionary runs are parametrized uniquely by the set of initial 
conditions given by its respective pair (Wo, <^o)- In establish- 
ing the starting models we follow Lupton & Gunn (1987) and 
take advantage of an almost always converging iteration pro- 
cedure between evaluation of the potential, brought into the 
form of a multipole expansion, and the determination of the 
density, where different types of Gauss-integrations are per- 
formed and thereby slightly different numerical procedures 
are employed (Einsel & Spurzem 1994). 

For simplicity and scaling reasons we choose the grav- 
itational constant G, the initial cluster mass Mi and the 
initial core radius to be one, i.e. to be our system units. 
The unit of time is expressed as (Cohn 1979) 



to = 



ri. [GMf 1 
GM ■ 47rr iV 



(15) 



F = 47r(Gm)^lnA, where A'^ is the total number of parti- 
cles. In A is the coulomb logarithm and m is the mean mass 
of particles (equal to 1/N in the present single mass mod- 
els). The first factor on the right hand side is simply the 
dynamical time of the system and evaluates to one with the 
units given above, while the second factor scales out with 
factors standing in front of the diffusion coefficients simpli- 
fying their expressions appreciably, then. An expression for 
the initial half mass relaxation time is taken from Spitzer & 
Hart (1971) 



0.138 



32.5 



N ■ 



G - m In A 



(16) 
(17) 



in units of to (Eq. (p^), and we assume A = 0.4 • N . Due to 
mass loss the relaxation time may change by a factor of 2-4 
during the runs presented here. 



3 NUMERICAL METHOD 
3.1 Computational grid 

In order to obtain an appropriate grid representation, e.g. a 
sufficiently resolved core and rectangular grid, a transforma- 
tion of the basic variables was performed. Following Cohn 
(1979) we take 



X{E) = - In 



2(i>c - Eq- E 



(18) 



where Eq is the energy of a circular orbit at the core radius 
of the cluster. Thus, inside the core (roughly <j>c ^ E ^ Eq) 
the relation is nearly linear, while the proportionality X oc 
— In |_E| in the halo improves the spacing of the radii of cir- 
cular orbits with given energies in the direction of the tidal 
boundary. For any specified energy the ^-component of the 
angular momentum is normalized to the maximum angular 
momentum for that energy, which again is the angular mo- 
mentum of the corresponding circular orbit as a function of 
E: 

For each time step rcirc{E) and Jz,o{E) are determined from 
the evolving potential in the equatorial plane by a simple 
Newton- Raphson scheme. At the same time the central po- 
tential changes, so that the complete {X, y)-grid is adapted 
to the new situation. 



3.2 Fokker-Planck step 

A finite difference scheme is used to integrate the diffu- 
sion part of the problem and the discretized equation sys- 
tem is solved implicitly with the help of a sparse matrix 
method (Henyey 1959) borrowed from a gaseous model code 
(Spurzem 1994, 1996). A Chang-Cooper scheme is applied to 
the energy direction in order to improve conservation char- 
acteristics (Chang & Cooper 1970). 

A double-logarithmic coordinate space grid is con- 
structed for the meridional plane {p,z), which extends ini- 
tially from 10"'^ core radii to somewhat beyond the tidal 
radius. During all runs described here it was adjusted to the 
contracting core radius rc{t) by the condition pi^Q — Zi=o = 
0.001 • rc{t) for the innermost grid point. This was usually 
done only about 4-6 times in each run in order to avoid 
rounding errors due to the necessary interpolation proce- 
dure involved. 

Because of deviations of the local velocity distribution 
from a rotating Maxwellian (Eq. (p^) especially in the outer 
halo, energy and angular momentum are not sufficiently con- 
served if the local foreground values of angular velocity and 
velocity dispersion are provided as parameters for the back- 
ground distribution function. Instead, these were used as 
starting values for an iteration procedure which determines 
new, consistent background values for lj and a under the 
condition that integrals of ft times the local first order dif- 
fusion coefficients < AE > and < AJz > (given in the 
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Appendix) over velocity space are zero (Goodman 1983), 
expressing the fact that locally no energy and angular mo- 
mentum may be generated or removed just by diffusion. 



3.3 Vlasov step 

After diffusion has taken place for an appropriate amount 
of time, the potential has to be advanced. The density is 
calculated from the distribution function according to 



i{p,z) = 



2tv 

P 



f{E,J,)dEdJ, 



(20) 



Given this density distribution, Eq. (Q) can be solved to ob- 
tain the potential. A multipole expansion is carried out to 
determine Dirichlet boundary values of (j) beyond the tidal 
radius in the meridional plane. Symmetry is assumed about 
the equatorial plane so that v. Neumann boundary condi- 
tions 



dz 



= 0' if 

z=0 Op 



= 



(21) 



p=0 



are applied at the inner boundaries. The system of equations 
is solved with the same sparse matrix method (Henyey et al. 
1959) as for the diffusion step, thereby taking account of the 
even sparser matrices occurring due to the Poisson equation. 

Having built the potential a new (_E, Jz )-grid is cre- 
ated, the meridional plane is scanned to determine the area 
A{E, Jz) and Eq. (^ is numerically integrated to obtain the 
adiabatic invariant q. It is required in order to fulfil the con- 
servation of f{q, Jz). At first, the new {E, Jz)-grid is trans- 
formed uniquely into a (q, Jz)-gnd. A second order bivariate 
Taylor expansion is then performed for each mesh point in 
the new grid to derive approximate values from the known 
values of f{q',J'z) in the old grid. I.e., f{E,Jz) is allowed 
to ch ang e. This again gives a new density distribution p in 
Eq. (EOl) and the procedure is repeated until convergence 
is reached. In doing this the firstly determined f{q,Jz) is 
retained as the old grid in order to save accuracy. Neverthe- 
less, this second order interpolation step turns out to be the 
relatively most important source of numerical error in the 
code. 

Goodman (1983) reports a test of the influence of the 
third integral on the procedure just stated. He relaxed an 
initial King model distribution function (non-rotating) vi- 
olently in assigning to it a Plummer potential. The result 
was a strong dynamical collapse with an increase in pc of a 
factor of 56. The collapse generated a flattening e — 0.055, 
which is attributed to the third integral, because there is 
no preferential axis of rotation or anisotropy initially, that 
could flatten the system. On the other hand, the test is very 
crude, because as Goodman himself states the above proce- 
dure is just first order in energy conservation between initial 
and final values of (ji and the collapse picture is not appro- 
priate to be described by adiabatic invariants. Though the 
effect of diffusion is to wipe out the (here artificially) gen- 
erated azimuthal anisotropy, we consider the flattening of 
nonrotating, isotropic conflgurations during their evolution 
as a more appropriate measure of the influence of the third 
integral. 



3.4 Computation and conservation characteristics 

The evolution of the system is followed numerically with 
time steps At chosen to be proportional to the central re- 
laxation time tr^ starting with At = 0.125 tr^ e.g. in case 
of a model with Wo = 6, but the coefficient was increased 
from time to time by a factor of | in order to have a frac- 
tional increase of the central density between 3—6%. In anal- 
ogy with the computations done by Cohn (1979) these runs 
ended with coefficients of about 2. Inbetween one Vlasov 
step, i.e. one recomputation of the potential, four Fokker- 
Planck steps, i.e. diffusion steps, were carried out. 

In order to stay as close to the King-model distribu- 
tion functions as possible we set up a strict tidal cutoff. The 
tidal boundary is adjusted so as to ensure that the mean 
cluster density is conserved throughout the evolution (see, 
e.g., Spitzer 1987). The fluxes at the Y = ±l-boundaries are 
set to zero as well as at the E = <^c-boundary. The only open 
boundary then is the E — _Etid-boundary, where the condi- 
tion /(-Etid, Jz) = was fixed. The first derivative of / with 
respect to E is non-zero at the boundary and is evaluated 
just inside the boundary in order to obtain accurate escape 
fluxes. With grid sizes of A^x = 100, Ny = 61, Np = Nz = 80 
used here the errors in mass, energy and angular momentum 
all accumulated to about 0.4%, 0.7% and 1.7% respectively 
for a typical model by the time the central density increased 
by about 5-6 orders of magnitude, when the runs were 
stopped. Goodman obtained errors 7.9% (12.7%) in mass, 
6.5% (2.9%) in energy and 16.5% (24.1%) in angular momen- 
tum for runs with Nx = 40, Ny = 10 {Nx = 20, Ny = 10) 
extending over 2i orders of magnitude, while Cohn (1979) 
reports an error in mass of only 0.04%, which may easily be 
explained by the application of spherical symmetry to real 
space in his simulations. 

Typical runs for the evolution of one model up to a den- 
sity increase of 5 orders of magnitude needed about 60 hours 
on a 100 Mhz hyperSPARC workstation (Maths Dept. in Ed- 
inburgh University) or 150 hours on a 40 Mhz SPARCstation 
(University of Kiel), respectively. Using a parallelized ver- 
sion of the code (message passing) on a CRAY T3D 40 hours 
were needed, when 8 processors were utilized. Although this 
speedup is not yet convincing, parallel processing is highly 
recommended for multi-mass versions of the present code. 

Although simulations with several choices of Wo have 
been carried out, the results presented here concentrate on 
those runs with Wo ~ 6 leaving the results with different 
Wo to the summary section and/or subsequent papers. Ad- 
ditional features indicated by those do not restrict or alter 
the conclusions drawn by using just the results from Wo = 6. 
The initial conditions of all models with Wo = 6 are summa- 
rized in Table |l]. Given are the rotation parameter uo, the 
ratio of rotational to total kinetic energy, the dynamical el- 
lipticity (c/. Eq. (p6[)), the ratios of tidal and half-mass radii 
to initial core radii, the central and finally the half-mass 
relaxation time in system units. 



4 NUMERICAL RESULTS 

4.1 Evolution of cluster structure 

The evolution of Lagrangian radii with time for models with 
Wo = 6.0 and luo ~ 0.0, 0.5, 1.0 is shown in Figure 0. In 
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Table 1. Initial conditions of all models presented with Wo = 6. 







edyn(O) 


rtid/rc{Oj 


rh/rdO) 


Trc(O) 


-f"rh(0) 


n no 


n on 
u.uu 


-U.UUl 


1 79 

-LO. I Z 


9 70 
Z. ( u 


19 24 


y 1 .OO 


0.05 


0.23 


0.002 


18.61 


2.70 


19.23 


91.77 




u.oy 


U.UIO 




z.uo 


19 22 


yu.ou 


n on 


Q QQ 

o.oo 


U.UOi 


10. OO 


9 


1 o on 


oy. / i 


0.30 


7.00 


0.105 


14.99 


2.62 


19.21 


87.73 


0.40 


11.23 


0.165 


13.08 


2.55 


19.22 


84.12 


0.50 


15.61 


0.224 


11.46 


2.48 


19.27 


80.49 


0.60 


19.81 


0.278 


9.94 


2.39 


19.40 


76.32 


0.70 


23.71 


0.327 


8.77 


2.30 


19.50 


71.78 


0.80 


27.18 


0.368 


7.69 


2.20 


19.71 


67.37 


0.90 


30.25 


0.403 


6.88 


2.12 


19.86 


63.24 


1.00 


32.99 


0.433 


6.22 


2.04 


20.02 


59.63 




2 



10 12 14 



tlx,,, 



Figure 1. Evolution of mass shells (Lagrange radii rss) for 
the model {Wo,i^o) = (6.0,0.60). Shown are the radii for mass 
columns containing the indicated percentage of total mass in the 
direction of the 6 = 54. 74° -angle. 
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Figure 2. Comparison of central density evolution for models 
with same Wo = 6.0, but different initial angular velocity param- 
eter uiQ = 0.0. ..1.0. The time is given in units of initial half mass 
relaxation times tr, ,. 



flattened systems the definition of Lagrangian radii is not 
straightforward. Just only for that purpose we assume here, 
since flattening is small, that deviations from spherical sym- 
metry are only up to second order in a Legendre expan- 
sion. Thus, we evaluate the Lagrange radii at a zenithal 
angle, where the effects of a probable flattening on the mass 
columns are expected to be less important: P2{cos6) — 0, 
which gives 6 = 54.74°. Thereby, we determine p(r,6) and 
then compute 

M(r) = I 47rp(r, 54.74° )r^dr . (22) 
Jo 

The tidal radius is determined from the condition 4'{p, z) = 
i?tid and is obtained here for the same angle as described 
above. 

The overall structure of evolution of mass shells closely 
follows what was given by previous non-rotating, isotropic 
or anisotropic models (e.g. Aarseth, Henon & Wielen 1974, 



Cohn 1979, Giersz & Heggie 1994a, b, Giersz & Spurzem 
1994) disregarding the existence of and evaporation through 
the tidal boundary. An expansion of mass shells larger than 
50% is observed and an initially smooth - later accelerated 
- contraction of the inner mass shells can be seen. The core 
radius decreases appreciably as well as the core mass. While 
no special sign of a gravo-gyro contraction (i.e. contraction, 
levelling off and (gravothermal) contraction again) is de- 
tectable, strong mass loss is evident in this diagram through 
the loss of complete mass shells. The tidal radius, which 
is not shown directly in the diagram, decreases in order to 
maintain the condition that the mean density in the cluster 
orbiting around a parent galaxy is conserved. 

Figure ^ shows a comparison of the density evolution of 
a sequence of models with the same initial Wo = 6.0 but with 
different wo = 0.0. ..1.0, respectively. The non-rotating model 
reaches singularity in its core after about tec ~ ILSltr^.o- 
This has to be compared with previously derived collapse 
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Figure 3. Comparison of scaled escape energy evolution for mod- 
els with same Wo = 6.0, but different initial angular velocity pa- 
rameter ujo = 0.0. ..1.0. The time is again given in units of initial 
half mass relaxation times tr^,i. 



time scales, which were in many cases determined for iso- 
lated Plummer models, giving tec ~ 15.5...17.6 tr^.o (Cohn 
1979, 1980, Takahashi 1996ab). Chernoff & Weinberg (1990) 
find times to reach core collapse for isolated King mod- 
els (without stellar evolution) with initial value Wo — 3.0 



{Wo = 7.0, Wo = 9.0) of tc 



{tc 



10.lt 



rh.O, 



tec = 2.23trh,o, respectively). Note that the latter two val- 
ues were not published in their paper, but reported by Quin- 
lan (1996), who himself finds a collapse time for the corre- 
sponding isotropic {f{E)) King model (case of Wo ~ 6.0) 
with tidal mass loss to be tec = 12.9trh,o. Detailed compar- 
ison with Quinlan's sequence of King models from Wo = 1 
up to Wo = 12 reveals a systematic difference of about 10 % 
between our and his results. An increase of the numerical 
resolution, i.e. the grid sizes, in the present models results 
in slightly larger collapse times thereby approaching Quin- 
lan's accurate values. Regarding the partly large differences 
between collapse times determined with distinct methods for 
(nearly) the same problems (isotropic or anisotropic mod- 
els) reported in the literature, we conclude that our results 
for the non-rotating models may be considered as to be in 
close agreement with previous work on that field. 

On the other hand Fig. |^ gives accurate informa- 
tion about the influence of rotation on collapse time. The 
strongest rotating model with tJ = 1.0 gives the smallest 
tec ~ 4.3trh,o, implying that rotation accelerates the col- 
lapse. While there is no three- or four-phase structure of 
contraction or collapse phases of different origin (Akiyama 
& Sugimoto 1989), the steeper slope in central density in 
the higher rotation case indicates the existence of a gravo- 
gyro contraction, which then transforms into gravothermal 
collapse which is more advanced than in the non-rotating 
case. 
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Figure 4. Evolution of central velocity dispersion versus cen- 
tral density for the same models as in Figure 3. The curve of 
the stronger rotating models show a remarkable upturn at the 
beginning (i.e. at low densities). 



A similar picture derives from Figure where the time 
evolution of the scaled escape energy 



xo 



(23) 



(cTc is the one dimensional central velocity dispersion) is 
shown. It is noteworthy that all models show the steep up- 
turn in xo, when they approach the collapse time and they 
arrive at values of xo between 8.0 and 9.0. This feature (and 
with this range of values) has already been described the- 
oretically by Lynden-Bell & Wood (1968) and numerically 
by Cohn (1979), so that we are in agreement with previous 
results interpreting this in terms of gravothermal instability. 
The effect of rotation can be seen in the earlier evolution: the 
model with oj = 1.0 is prevented from the usual deepening 
of the potential scaled by the velocity dispersion in the cen- 
tre of the cluster and stagnates near its starting value of xo- 
Considering the stronger increase in central density inferred 
from Figure ^ for this model at the same time, the reason for 
the behaviour in Fig. ^ will be found in a relatively strong 
increase in the central velocity dispersion. The sources of 
this effect then subside slowly and a smooth turnover into 
the fast collapse is traced, and so this model crosses the 
paths of the still gently contracting slow- and non-rotating 
conflgurations. 
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Table 2. collapse parameters of all simulated models with Wo • 
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Figure 5. Evolution of angular momentum versus central density 
for all rotating models with = 6.0. 



To clarify these processes we plot in Figure ^ the cen- 
tral velocity dispersion against central density. The evolu- 
tion starts in the lower left corner of the diagram. While 
the slowly and non-rotating models show the typical sign 
of an initial reluctance to increase in velocity dispersion as 
is the case for isotropic and anisotropic systems (e.g. Cohn 
1979), the rapidly rotating model {uio = 1.0) seems to suf- 
fer a significant thermalization and heating of its core. The 
evolutionary paths then join again and start the self-similar 
collapse phase, which is characterized by the quantity 



7 = 



dlntJc 
din pc 



(24) 



which represents twice the slope in the double logarithmic 
plot ((Jc, not (Tc is plotted). From all models of our 2D- 
calculations we derive an average 7 = 0.109 with no de- 
pendence of 7 on the amount of initial rotation, implying 
that initial rotation seems to have a minor influence on col- 
lapse characteristics. Cohn's (1979) and Takahashi's (1996a) 
2D-results show 7 — 0.12 and 0.10, respectively, and Cohn 
(1980) obtains 7 = 0.10 from his isotropic systems in rough 
agreement with our result. Goodman (1983) also derives a 
value of 0.10 for 7 with his low resolution grids. 

The gravo-gyro picture by Hachisu (1979) indicates that 
the positive variation of temperature in the core is due to 
a rise in entropy caused by viscosity effects, when angu- 
lar momentum is transferred to disordered motion, and is 
thereby effectively carried outwards (negative specific mo- 
ment of inertia). In Fig. |^, a plot of central angular velocity 
versus central density is given in the same manner as the 
log((Tc) — log(pc) diagram above, but only for the rotating 
models. There is a pronounced initial increase in angular 
velocity for the most rapidly rotating model, which after 
about 2.5 trf^,i reaches a plateau and finally turns over into 
a power law relation of constant slope in this logarithmically 
plotted diagram when roughly 4 fr^.i have passed by. The 
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increase in loc is connected with the heating process evident 
from Figure ^ 

The stagnation phase after the initial increase marks 
the time, when most of the total angular momentum is di- 
minished (transported into the outer halo or evaporated) 
ending the state of negative specific moment of inertia in 
the core. The further evolution is now determined by self- 
similar evolution in the gravothermal picture, as can be seen 
if we take a rough measure of the quantity 

S^^. (25) 
din pc 

Inspection of Figure |5| reveals 5 = 0.06... 0.08 (with a small 
dependence on the amount of initial rotation), which means, 
that azimuthal (rotational) velocities in the decreasing core 
increase in roughly the same way as any other velocity (in- 
ferred from 7/2 ~ 0.055). Goodman found values 5 = 0.14 
for grid size {Nx x Ny) = (20 x 10) and 5 = 0.10 for a 
(40 X 10)-grid. 

Several useful quantities derived from the models pre- 
sented here are shown in Tab. ^. In column (1) the initial 
angular velocity in system units is given {Wo = 6), col. (2) 
gives the collapse time, col. (3) the collapse rate, col. (4) 
the number of current central relaxation times, Trcm, until 
complete collapse in the self-similar evolution phase, col. (5) 
and (6) the exponents used in the corresponding equations of 
state of the core, while col. (7) and (8) state the percentage 
of initial mass and angular momentum loss per half mass re- 
laxation time. Note, that the non-rotating model gives mass 
loss rates in full agreement with those values reported in the 
literature (e.g. Spitzer 1987), with strongly increasing val- 
ues towards models of higher rotation. Unfortunately, the 
well known value of 320 for Trcm (Cohn 1980) is underesti- 
mated by our models, which indicates that these suffer from 
worsening numerical resolution, when the collapse proceeds 
further and further. 



© 1997 RAS, MNRAS 000, 



Dynamics of rotating stellar systems 9 




2.0 4.0 6.0 



10.0 12.0 14.( 



tlx., 



Figure 6. Evolution of total mass retained in the cluster models 
with Wo = 6.0. 



4.2 Evolution of cluster shapes 

As has already been shown in Figure 0, our models - subject 
to the strict tidal boundary formulation given above - suffer 
strong mass loss, which can also be seen from Figure |^. Here, 
a comparison of preserved mass as a function of time be- 
tween the models with same Wo = 6.0 but different rotation 
{cuo.i ~ 0.0. ..1.0) is given. In general this diagram reveals 
information about evaporation of angular momentum from 
the halo, which drives cluster evolution with respect to rota- 
tion. As expected the fastest rotating model shows strongest 
mass loss indicated by the slope of the curves, and within 
its evolution an even steeper mass loss can be found in the 
initial two half mass relaxation times, which were associated 
with the gravogyro phase of the system. 

By transforming formulae derived by Binney (1978), 
who applies the tensor- virial theorem to a system composed 
of similarly situated ellipsoids, Goodman (1983) defines the 
so-called dynamical ellipticity edyn, which then is given by 



2Trot + STo-^ — To- (1 + 2s ) arccos s — 3s\/l — 



arccos s 



(26) 



where s = b/a = 1 — edyn is the axis ratio of the ellipsoids, 
which are assumed to be oblate spheroids (i.e. not triaxial). 
Trot is the rotational energy, To-^ is the energy contained in 
the azimuthal component of the velocity dispersion and To- 
is the energy associated with all components of the velocity 
dispersion. For sufficiently small edyn this relation may be 
expressed as 



8 

J. edyn . 



(27) 



The evolution with time of edyn (calculated from Eq. (p6|)) 
is plotted in Figure |^ for the three models with Wo = 6.0. A 
steep decrease in ellipticity can be seen for the initially most 
strongly rotating model and the final states of all models 



lack significant fiattening. Due to less effective mass loss, the 
angular momentum transport beyond the tidal boundary is 
smaller in the more moderately rotating models, so that the 
respective curves are able to cross each other. 

Spherical symmetry is not exactly preserved for the 
non-rotating model (see Fig. 7, curve uio =0), which can 
be explained by the creation of anisotropy in its halo, when 
azimuthal pressure is compared with that within the merid- 
ional plane. Because there is a preferential direction due 
to the chosen coordinate system, i. e. the z-direction, no 
anisotropy may develop in planes tangential to the poles of 
the system. Thus, just as radial anisotropy develops in gen- 
eral spherical systems (e.g., Cohn 1979, Louis & Spurzem 
1991), it evolves in the equatorial plane of our model con- 
figurations, but not along the z-axis with p <^ rc, which 
is caused by the neglect of the third integral. Moreover 
the present formulation of the tidal boundary produces an 
anisotropy profile, which shows the usual rise towards the 
halo starting from the isotropic core, but then, as it ap- 
proaches the tidal boundary, falls below zero, because stars 
on radial orbits are evaporated more effectively than stars 
on tangential orbits: the former easily gain energy in core 
passages, so that those orbits are deplenished. On the other 
side stars visiting regions near the poles are all coupled to 
the core, because no third integral specified prevents them 
from doing that (e.g., in the spherical symmetric case, if 
it is non-zero for the orbit given). 

Neglect of the third integral also has effects on orbits 
when slow potential changes take place. Consider circular 
orbits in the meridional plane and equatorial plane and as- 
sume a currently spherical potential. A gentle contraction of 
the system will retain a circular orbit in the equatorial plane, 
since it is represented accurately by its [E, J2)-pair (it will 
just be shifted to that energy, where Jz,o{E) = J°''' for 
angular momentum has to be conserved), while the former 
orbit is not distinguished by our coordinate systems from 
radial orbits of the same energy. Therefore, all orbits with 
Jz ~ are shifted to new energies (via the condition of adia- 
batic invariant conservation) corresponding to the ensemble 
mean with respect to third integral, and this may not be the 
same energy in the meridional circular orbit case as that of 
the equatorial circular orbit. In particular, the probability 
of retaining a circular orbit is not necessarily one. 

Nonetheless, inspection of Figure |^ indicates, that the 
ellipticity of the non-rotating model rises due to the neglect 
of the third integral only to values edyn ~ 0.01, well below 
the uncertainties inferred from observations. 

In Figure |^ the variation of unprojected ellipticities 
e = 1— 6/a with radius is shown for the model with ljo ~ 0.50 
for several evolutionary times. The current core radii are in- 
dicated by triangles on the curves, respectively. The elliptic- 
ities of the core are evidently smaller as is the case already 
for the starting models. The ellipticity peaks at about the 
half mass radius where it is appreciably larger than in the 
core but with a nearly constant offset (Ae « 0.07 in this 
case). This indicates, that we indeed expect a significant el- 
lipticity variation within globular clusters as has been stated 
by Geyer et al. (1983). It is remarkable that the contraction 
or decrease of the core is not refiected in this picture. The 
horizontal lines in this Figure give for comparison the re- 
spective dynamical ellipticities edyn for that time step. It 
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Figure 7. Evolution of dynamical ellipticity edyn as defined by 
Goodman (1983) (Eq. (27) see text for explanation) for all models 
with Wo = 6.0. Fig. a) gives the evolution versus time measured 
in units of the half mass relaxation time, while Fig. b) shows the 
evolution versus the scaled escape energy xq. 



seems, that this quantity represents very well a global mean 
unprojected ellipticity. 

Strong evidence for rotation in globular clusters had 
been derived from the observed coincidence of ellipticity 
and rotational velocity profiles, both being scaled arbitrarily 
(Meylan & Mayor, 1986). In Figure ^ we show for compari- 
son for the same model and time steps as in Fig. ^ the corre- 
sponding profiles of rotational velocity. Prior to gravother- 
mal instability a coincidence of the profiles just stated may 
indeed be watched, while for the rotational velocity curve 
corresponding to the gravothermal collapse phase 

Recently, maps of rotational velocities and azimuthal 
as well as meridional velocity dispersions have been con- 
structed from non-parametric fitting of large samples of 
kinematical data for uj Cen (Merritt et al. 1997). Disre- 
garding for the moment the simplicity of the present models 
(i.e. single mass models, no stellar evolution, etc.) we chose 
a model from our simulations, which roughly reproduces the 
ellipticity and concentration parameter (0.12 and 1.36, re- 
spectively) derived for uo Cen from observations. The veloc- 
ity maps obtained (Wo = 6, uo = 0.6, t = 3.77 rrhg) are 
shown in Fig. |l^. These maps agree with those from Merritt 
et al. (1997) in the morphological structure, e.g. the oblate 
isovelocity contours in the case of or the torus-like con- 
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Figure 8. profile of unprojected ellipticity for model with Wq = 
6.0 and uiq = 0.50 for several evolutionary times as indicated. The 
pairs in brackets give xq and t/t^i^, respectively. The horizontal 
lines on the right side indicate the corresponding dynamical el- 
lipticities. 
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Figure 9. profile of rotational velocities v^ot for model with Wq = 
6.0 and uiq = 0.50 for several evolutionary times as indicated. The 
pairs in brackets give xg and t/ty-i^, respectively. The sequence of 
curves starts with the pair (0.00, 6.21) at the top. The actual core 
radii are denoted by the triangles put on the curves, respectively. 
The core radius of the last model plotted is situated left from the 
figure area. 
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Figure 10. Contour maps a) of velocity dispersion in the merid- 
ional plane b) of the velocity dispersion in azimuthal direction 
and c) the rotational velocity for a model with Wo = 6, loq = 0.6 
and t = Trhn ■ 



tours when the rotational velocity map is considered. Note, 
that the axes are measured in units of initial core radii and 
that the current core radius has already decreased to about 
0.44 re-. 



5 SUMMARY AND DISCUSSION 

We have performed 2D-Fokker-Planck simulations modelling 
the evolution of rotating stellar systems. The main results 
can be summarized as follows: large amounts of initial rota- 
tion drive the system into a phase of strong mass loss while 
it contracts slightly. The core is rotating even faster than be- 
fore although angular momentum is transported outwards. 
At the same time the core is heating. Given these features 
we associate this phase with the gravo-gyro phase found 
by Hachisu (1979). The total collapse time is shortened ap- 
preciably by this effect, but the increase in central angular 
momentum levels off after about 2 — 3 fr^.i indicating that 
the source of this 'catastrophe' ceases, i.e. it is not really 
a catastophe. Finally, the central angular velocity increases 
again, but with a rather small power of the central density - 
nearly the same power as for the central velocity dispersion 
during self-similar contraction. We propose to search for a 
self-similar solution of rotating cluster models in the future, 
which should agree with our findings. 

Due to the non-locality of relaxation, the processes of 
mass loss, angular momentum transport and gravo-gyro and 
gravothermal 'catastrophe' in our models are not easily dis- 
entangled and as well are not entirely the same as e.g. in ro- 
tating gas spheroids. Additionally, even in that case no the- 
oretical formulation of the problem has been found, yet. The 
only theoretical investigations undertaken so far concern ro- 
tating gas cylinders (Inagaki & Hachisu 1978, Hachisu 1979) 
and rotating gas disks (Hachisu 1982). But the mechanisms 
taking place in our comparatively complicated structure of 
a slightly flattened spheroid may be identified by analogy 
with an idealized transition between the cylinder and disk 
configurations of Hachisu (1982). The mechanism of gravo- 
gyro 'catastrophe' had been introduced by Hachisu (1979) 
for rotating gas cylinders, but these models do not comprise 
gravothermal 'catastrophe' for reasons concerning the cur- 
vature of the metric configuration (Hachisu 1982). On the 
other hand the investigation of interactions between grav- 
ogyro and gravothermal processes may be accomplished by 
considering the disk configuration in which both processes 
are active. There, it is shown that those terms in the tensor, 
which describes the hydrostatic readjustment of the system 
to a variation, that denote a coupling of quantities (e.g. an- 
gular momentum with temperature) play an important role 
for further evolution, thereby, e.g., adjusting the gravogyro 
to the gravothermal contraction. Then, the time scale of evo- 
lution will be given just by the efficiency of heat transport. 

The time needed to reach collapse for each model is 
consistent with the individual mass loss rate (see Fig. (??), 
which itself is connected with the amount of angular momen- 
tum transport due to viscosity. The larger the initial amount 
of rotation the stronger the mass loss (see Fig. <E>??). Con- 
trary to the simple models of Agekian (1958) stellar esca- 
pers from the system reach their escape energy in encounters 
mostly inside or near the core. 

Due to both features - the total time to reach collapse 
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depends on the amount of mass loss induced by viscosity 
effects and the solution for the gravothermal collapse is only 
moderately modified by gravogyro effects - it may be con- 
cluded, that the gravogyro 'catastrophe' is coupled to the 
gravothermal contraction and therefore does not evolve on 
time scales of angular momentum transport (viscosity co- 
efficient) but rather on time scales of heat transport. This 
is clearly indicated by the lack of an early, separate central 
contraction phase of Lagrangian radii (see Fig. (??)) as it 
seemed to be implied by the N-body simulations of Akiyama 
& Sugimoto (1989). Note, however, that they did not include 
a tidal boundary. On the other hand gravogyro effects influ- 
ence cluster evolution indirectly by increasing the mass loss 
activity thereby reducing the total collapse time. 

One of the main results of this work is that a cluster 
system presently comprised of nearly spherical clusters may 
have had a distribution of cluster ellipticities up to strongly 
flattened, rotating stellar systems. This raises question for 
the origin of globular clusters: the picture presented here 
would be fully consistent with a formation from (obviously 
significantly rotating) giant molecular cloud like structures, 
when comparing their mass and angular momentum content 
(Akiyama & Sugimoto 1989). Other mechanisms, which are 
still not in contradiction to our results, incorporate the for- 
mation of rotating globular clusters from previously exist- 
ing binary clusters, which suffer a synchronization instabil- 
ity leading to a merger (Sugimoto & Makino 1989). It is 
found that employing this mechanism a maximum fiatten- 
ing e ~ 0.3 of the globular clusters can be explained. The 
advantage of this model is that it explains very well several 
features of the Magellanic Clouds globular cluster system. 
The stronger flattening of Magellanic Cloud clusters as com- 
pared to those of the galaxy is in the light of our simulations 
understandable as a consequence of their dynamical age. An- 
other possible mechanism to form binary clusters is collapse 
from a shell (Theis 1996) . The formation of (single) clusters 
from collapsing shells has been proposed by Brown et al. 
(1991, 1995). 

The rotation curves decrease with time, as do the ellip- 
ticity profiles, but the relative shapes of the latter do not 
vary much. The maximum values of rotational velocity and 
ellipticity remain at about the half mass radius throughout 
the evolution. In accord with our results Gebhardt et al. 
(1995) observe in 47Tuc by using non-parametric techniques 
an increase in the rotational velocity towards the half-mass 
radius and nearly solid body rotation inside the core. How- 
ever, using the rotational velocity derived from integrated 
light (Gebhardt et al. 1994) they obtain an increasing angu- 
lar velocity for the region inside half a core radius. Moreover, 
47 Tuc reveals an already advanced stage of evolution, so 
that the results presented here indicate that a maximum in 
the velocity curve occurs at tens of current core radii (pre- 
liminary modelling gave about 20rc(t) for a 47 Tuc - like 
model cluster). While the first argument may be reconciled 
with the results presented here (solid body rotation in the 
core and a rotation velocity peak at about the half-mass ra- 
dius, which is roughly 10 core radii in the case of 47Tuc), 
when the bad statistics of the data beyond ~ 2rc in the 
observations are considered, the rise in angular velocity in 
the innermost central parts of the core, which are indicated 
by those observations may not be explained by our actual 
single mass models. There are several possible reasons, why 



our models may not yet be appropiate for 47Tuc: i) in 47Tuc 
as in any other cluster near core bounce there is significant 
activity of hard binaries going on, whose reaction products 
have been observed (Meylan et al. 1991); the present mod- 
els do not yet incorporate binary energy generation; ii) the 
cluster mass function is observed to vary with radius and 
time both in spherical models and observations. Models in- 
cluding mass segregation effects may differ significantly from 
the single-mass case, as it is known for non-rotating clus- 
ters. For example it may be speculated that in deep core 
collapse high mass stars, which quickly segregate towards 
the centre, cannot loose angular momentum as efficiently as 
the stars in our simple models with equal masses do. iii) 
Finally, due to its very peculiar structure 47Tuc could be 
a non-axisymmetric cluster due to a recent encounter with 
the galactic disk, bulge or other cluster. 

Future work, which will be published in the next pa- 
pers of this series, consists of the incorporation of an energy 
source due to formation and hardening of three-body bi- 
naries (Bettwieser & Sugimoto 1984, Heggie & Ramamani 
1989, Lee, Fahlman & Richer 1991) in order to proceed into 
the post-collapse phase, and the yet completely uninvesti- 
gated influence of differently rotating mass groups. Mass 
segregation would presumably alter the efficiency of angular 
momentum transport significantly. Such models will allow 
more detailed modelling of observed clusters, such as 47 Tuc 
or others. 

Also it is important to consider the effects of different 
tidal boundary conditions and the dynamical infiuence of 
mass loss due to stellar evolution, both of which can al- 
ter the cluster structure considerably. As for the treatment 
of the tidal boundary it affects e.g. the anisotropy profile 
near the boundary. Drukier's (1995) boundary formulation 
with delayed mass loss beyond the tidal boundary, respec- 
tively, are worth trying out in this context. Finally, since the 
gravogyro contraction seems to take place in the very early 
evolution of the cluster, different initial conditions also may 
have important bearings on the rotational structure of the 
system. 

We stress the necessity to perform N-body simulations 
of rotating systems in order to indicate the correctness of 
the assumptions made in the present Fokker-Planck models, 
especially to see the effects due to the neglect of the third 
integral. 
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APPENDIX A: DERIVATION OF THE FLUX 
COEFFICIENTS 

When deriving diffusion coefficients several ways may be 
chosen for doing it. We decided to follow the approach of 
Rosenbluth, MacDonald & Judd (1957) involving covariant 
derivatives of tensorial objects. At first, the coordinate sys- 
tems will be introduced. While we have cylindrical coordi- 
nates in coordinate space {p,z,ip), the following symmetry 
is applied to velocity space: = v = {Vp + v'^ + vl)'^; 

= ip = arctan(«p/ii2); q^ — v^, where (vp,Vy,,Vz) are lo- 
cal cartesian velocity coordinates. The corresponding metric 
tensor (a'"") then reads as 



1 





V 





(-^+"^) 





V 





1 



(Al) 



with the volume element a := det(afj^) = 1/ det(a'"^) = 
v'^. The tensorial form of the Fokker-Planck equation may 
generally be written as 

f:^ = -ifT!:\. +i(,fsr),M. , (A2) 

where the commas denote covariant derivatives, and the 
subscript a indicates particle species. The factor Ta ~ 
ATvGrria In A contains the usual Coulomb logarythm. The dif- 
fusion coefficients of Cartesian coordinate systems may then 
be expressed as tensorial objects 

^ * " a^^iha),. (A3) 



< Av'' >a= 



(A4) 



The functions h and g are the so called Rosenbluth poten- 
tials: 



ha{v) 



nig + nib 
mt 



dvfMvf)- 



g{v)^T,b / difffb{vf)\v - iff\ 



(A5) 



(A6) 



After some lengthy calculations involving Christoffel 
symbols we arrive at the following expressions for the tensors 
given above (note, that symmetry is assumed about ?/)): 



n. 



dh 



V dv, 



dh 

dv 





dh ^ dh 

V dv dvu> 
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dv^ V dvdvtr 
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n3 Qy 

Thus, employing the relations 



+ 2- 



(AlO) 
(All) 
(A12) 

(A13) 
(A14) 
(A15) 

(A16) 
(A17) 



with the r symbol denoting a Christoffel symbol of the sec- 
ond kind, here, the Fokker-Planck equation consists of terms 



and 
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dv 
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~^ V ] dv 
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< AE >= E,^ < Av" > +-E,^, < Av^Av" > (A22) 



< AJ, >= J,f, < Av" > +- J,M. < Av^Av'^ > (A23) 



< {AEf E,^ E,, < Av^Av" > 



< {AJ,y >= J,p J„ < Av^Av" > 



< AEAJ, >= E,^ J„ < Av'^Av'^ > 



(A24) 
(A25) 
(A26) 



If we now assume an isotropic background distribution 
f{vf), integrals h and g are easily simplified (cf. Spitzer 
1987, p. 36). For this purpose, we express the velocity v (in- 
ertial frame) in terms of velocities {u + pQe^) in the corotat- 
ing frame. Q is the angular velocity of the corotating frame. 
Thus, f{u) is isotropic and the derivatives of h and g with 
respect to v and v^, ^j^, etc.) must be transformed to 

those with respect to u only. For example, ~ 



Inserting eqs. ( [AT] ) to ( [AISD in eqs. (|A20|), (|A2l[) and em- 
ploying these results again in eqs. (A22...4.26) under con- 
sideration of the transformation to the corotating frame, we 
arrive at 



<AE> = 



u 



u 



dh 
du 



~^ 2dv? 
u du 



< AJ, > = — 



J, f/£l\ dh 
u u I du 



< (AEY > 



u + 2J-M - 2p^Q.^ + 



(A27) 



(A28) 



Most terms involving S vanisch due to further derivations 
with respect to the coordinate of symm etry, tp, but one of 
them is retained in the last term in Eq. (A19), therein writ- 



ten explicitly. The diffusion coefficients in transformations 
between curvilinear coordinate systems may now be identi- 
fied using 



1 



(A20) 
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< Av^'Av'' > = sr 

i a 



(A21) 



where the additional term in Eq. (A20) with respect to the 
case of E q. (A3) originates from the second (last) term of 
Eq. (|aT9|). 



It is convenient to treat the problem in energy - angular 
momentum - space such that the diffusion coefficients just 
derived have to be transformed to the new velocity variables 
E = + "l>(p, z) and = pv^p. This can be accomplished 
by using the following simple formula: 



< {AJ.f > 



du^ 



— + 
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< AEAJ, > 
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Dee = — {2uJzn-up'^n'^ +u^) Ei{v 



, r, I 2^2 Jz! 
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The Fokker-Planck equation is usually recast in flux conser- 
vation form: 



Df _ OFe dFj^ 



Dt dE dJ, 



(A32) 



where we still have not applied an orbit average. The term 
on the left hand side represents the Vlasov-term of the full 
coUisional Boltzmann equation. The fluxes F are given by 



+ — + - ]F4{u), 



u u 
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Fe = -Dee% - Dej. If - De/ 
oE oJz 



F,.=-D,^,.^^-D,J^^D,J. 



(A33) 
(A34) 



Comparing eqs. (A32 - A34 ) with the original form (e.g. 
Spitzer 1987) the flux coefflcients D may be identified to 

Dee = \< (AEf > 
Dej, = ^ < AEAJz > 



Dej, 



u u 
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Inserting the expressions for the diffusion coefficients given 
above, one obtains finally the desired flux coefficients for the 
axially symmetric application with rotation. We find 



De = 4nl^-^+u]F2{n), 

u U I <■ ' 



D, 



= 4n [ ^ - ^ ] F2{u), 
'it u ' 
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The functions Ei and Fi are constituent parts of the Rosen- 
bluth potentials and their derivatives. 



E, 
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u' 
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u"' f{u')du' 
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(A42) 



Given an expression for the background distribution func- 
tion (section |2|, rotating King-model used throughout this 
paper) these functions are evaluated for each position in the 
meridional plane, when the local density, mean particle ve- 
locity and mean particle angular velocity are specified. 

The next step is to orbit average these fiux coefficients 
and again to transform them to {X, F)-coordinates used 
throughout our numerical simulations. 
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